suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(ChIPseeker))
suppressPackageStartupMessages(library(clusterProfiler))
suppressPackageStartupMessages(library(magrittr))
suppressPackageStartupMessages(library(ggvenn))

small_legend <- function(.plot, pointSize = 1, textSize = 8, spaceLegend = 0.8) {
  .plot <- .plot +
    guides(shape = guide_legend(override.aes = list(size = pointSize)),
           color = guide_legend(override.aes = list(size = pointSize))) +
    theme(legend.title = element_text(size = textSize),
          legend.text  = element_text(size = textSize),
          legend.key.size = unit(spaceLegend, "lines"))
  return(.plot)
}

suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene))
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
annoDb <- "org.Hs.eg.db"

files <- c("CN"="files/CN.peaks.xls", "SH"="files/SH.peaks.xls")

chr_list <- read.table(files[1], header=TRUE) %>% distinct(chr) %>% pull(chr) %>% as.vector()
chr_list <- chr_list[nchar(chr_list) <= 5]


peak <- readPeakFile(files[1])
peak <- peak[seqnames(peak) %in% chr_list]

covplot(peak, weightCol="pileup")

peak <- readPeakFile(files[2])
peak <- peak[seqnames(peak) %in% chr_list]

covplot(peak, weightCol="pileup")

# covplot for all peak files
# ...

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)

tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), facet="row") + ggsci::scale_color_aaas()
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500, facet="row")  + ggsci::scale_color_aaas() + ggsci::scale_fill_aaas()


tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=ggsci::pal_aaas()(length(files)))


peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb, tssRegion=c(-3000, 3000), verbose=FALSE)

plotAnnoBar(peakAnnoList) %>%  small_legend()

plotDistToTSS(peakAnnoList) %>%  small_legend()



genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
ggvenn(genes, fill_color = pal_aaas()(length(files)))


compMF <- compareCluster(geneCluster   = genes,
                         fun           = "enrichGO",
                         ont           = "MF",
                         pvalueCutoff  = 0.05,
                         OrgDb         = annoDb,
                         pAdjustMethod = "BH")

dotplot(compMF, showCategory = 10, title = "Molecular Function Enrichment")

compBP <- compareCluster(geneCluster   = genes,
                         fun           = "enrichGO",
                         ont           = "BP",
                         pvalueCutoff  = 0.05,
                         OrgDb         = annoDb,
                         pAdjustMethod = "BH")

dotplot(compBP, showCategory = 10, title = "Biological Process Enrichment")

compCC <- compareCluster(geneCluster   = genes,
                         fun           = "enrichGO",
                         ont           = "CC",
                         pvalueCutoff  = 0.05,
                         OrgDb         = annoDb,
                         pAdjustMethod = "BH")

dotplot(compCC, showCategory = 10, title = "Cellular Component Enrichment")


compKEGG <- compareCluster(geneCluster   = genes,
                           organism      = "hsa",
                           fun           = "enrichKEGG",
                           pvalueCutoff  = 0.05,
                           pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 10, title = "KEGG Pathway Enrichment Analysis")